data589project_clean

DATA import

library(spatstat)

## Loading required package: spatstat.data

## Loading required package: spatstat.geom

## spatstat.geom 3.2-9

## Loading required package: spatstat.random

## spatstat.random 3.2-3

## Loading required package: spatstat.explore

## Loading required package: nlme

## spatstat.explore 3.2-7

## Loading required package: spatstat.model

## Loading required package: rpart

## spatstat.model 3.2-11

## Loading required package: spatstat.linnet

## spatstat.linnet 3.1-5

## 
## spatstat 3.0-8 
## For an introduction to spatstat, type 'beginner'

library(sp)
library(sf)

## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE

library(splines)
load("wolf_clean_nondup.RData")
load("BC_Covariates.Rda")

#First moment descriptive statistics

General distribution

plot(wolf_ppp, main='Coyotes')

## Quadracount

Q <- quadratcount(wolf_ppp, nx = 5, ny = 5)

plot(wolf_ppp,
     pch = 20,      
     cex = 0.5,    
     main = "Quadrats Visualization")

plot(Q, add = TRUE, col = "red")

## Quadratest

quadrat.test(Q)

## Warning: Some expected counts are small; chi^2 approximation may be inaccurate

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  
## X2 = 1319.3, df = 20, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 21 tiles (irregular windows)

Intensity

Q <- quadratcount(wolf_ppp, nx = 5, ny = 5)

plot(intensity(Q, image = TRUE), main = "Intensity", col = terrain.colors(256))

# Overlay the 'wolf_ppp' point pattern on the density plot
plot(wolf_ppp,
     pch = 20,        # Type of point, here '20' represents a filled circle.
     cex = 0.6,       # Point size.
     col = "blue",    # A single color for all points, you can choose your preferred color.
     add = TRUE)      # Indicates that points should be added to the existing plot.

## Heatmap

wolf_density <- density(wolf_ppp)

# Plot the density estimate
plot(wolf_density, main = "Heat map")

plot(wolf_ppp,
     pch = 20,        # Type of point, here '20' represents a filled circle.
     cex = 0.5,       # Point size.
     col = "green",    # Color of points, choose your preferred color.
     add = TRUE)      # Indicates that points should be added to the existing plot.

From above preliminary analysis, we find the distribution is inhomogeneous.

Covariates Analysis

Elevation <- BC_wolf$Elevation
Forest<-BC_wolf$Forest
HFI <- BC_wolf$HFI
Dist_Water <-BC_wolf$Dist_Water

We cut the 4 covariates into 5 classes and explore their relationship with coyote distribution.

Elevation

elev_cut <- cut(Elevation,
                5,
                labels = c("very low", "low","medium","high","very high"))
# plot the elevation class image
plot(elev_cut,
     main = "Elevation classes")
# overlay the park locations
plot(wolf_ppp,
     cols = "black",
     pch = 16,
     cex = 0.6,
     add = TRUE)
plot(wolf_ppp,
     cols = "white",
     pch = 16,
     cex = 0.5,
     add = TRUE)

wolf_elev_class <- elev_cut[wolf_ppp]
# Create a table of counts within each class
table(wolf_elev_class)

## wolf_elev_class
##  very low       low    medium      high very high 
##       233       108         3         0         0

Forest Cover

forest_cut <- cut(Forest,
                5,
                labels = c("very low", "low","medium","high","very high"))
# plot the elevation class image
plot(elev_cut,
     main = "Forest classes")
# overlay 
plot(wolf_ppp,
     cols = "black",
     pch = 16,
     cex = 0.6,
     add = TRUE)
plot(wolf_ppp,
     cols = "white",
     pch = 16,
     cex = 0.5,
     add = TRUE)

wolf_forest_class <- forest_cut[wolf_ppp]
# Create a table of counts within each class
table(wolf_forest_class)

## wolf_forest_class
##  very low       low    medium      high very high 
##       165        83        22        52        19

HFI

HFI_cut <- cut(HFI,
                5,
                labels = c("very low", "low","medium","high","very high"))
# plot the elevation class image
plot(HFI_cut,
     main = "HFI classes")
# overlay
plot(wolf_ppp,
     cols = "black",
     pch = 16,
     cex = 0.6,
     add = TRUE)
plot(wolf_ppp,
     cols = "white",
     pch = 16,
     cex = 0.5,
     add = TRUE)

wolf_HFI_class <- HFI_cut[wolf_ppp]
# Create a table of counts within each class
table(wolf_HFI_class)

## wolf_HFI_class
##  very low       low    medium      high very high 
##        19        38        73        61       127

Distance to water

Dist_Water_cut <- cut(Dist_Water,
                5,
                labels = c("very close", "close","average","far","very far"))
# plot the elevation class image
plot(Dist_Water_cut,
     main = "Dist_Water classes")
# overlay
plot(wolf_ppp,
     cols = "black",
     pch = 16,
     cex = 0.6,
     add = TRUE)
plot(wolf_ppp,
     cols = "white",
     pch = 16,
     cex = 0.5,
     add = TRUE)

Dist_Water_class <- Dist_Water_cut[wolf_ppp]
# Create a table of counts within each class
table(Dist_Water_class)

## Dist_Water_class
## very close      close    average        far   very far 
##        332         13          2          0          0

Rhohat estimate

rho_elev <- rhohat(wolf_ppp, Elevation)
rho_forest <- rhohat(wolf_ppp, Forest)
rho_hfi <- rhohat(wolf_ppp, HFI)
rho_water <- rhohat(wolf_ppp, Dist_Water)

plot(rho_elev, xlim=c(0, max(Elevation)), main="ρ vs Elevation")

Lower elevations have a higher than average density of coyote occurrences.

plot(rho_forest, xlim=c(0, max(Forest)), main="ρ vs Forest")

When forest cover is <50%, the density seems higher.

plot(rho_hfi, xlim=c(0, max(HFI)), main="ρ vs HFI")

The intensity of sightings is initially low and increases sharply as the HFI approaches 1.

plot(rho_water, xlim=c(0, max(Dist_Water)), main="ρ vs Dist_Water")

The relationship does not seem to be linear or simply exponential. The significant width of the confidence interval in the latter part of the plot for the distance from water indicates there is considerable uncertainty in the intensity estimates at greater distances.

From above, it seems all covariates have different relationship with coyote distribution.

Second moment descriptive statistics

Morisita’s index

W <- owin(xrange = c(min(wolf_clean_unique$X), max(wolf_clean_unique$X)), yrange = c(min(wolf_clean_unique$Y), max(wolf_clean_unique$Y)))

wolf_ppp_2 <- ppp(x = wolf_clean_unique$X , y = wolf_clean_unique$Y, window = W)

miplot(wolf_ppp_2, main = "", pch = 16, col = "#046C9A")

Zoom in

miplot(wolf_ppp_2, main = "", xlim = c(0, 60000), pch = 16, col = "#046C9A")

Under the assumption of homogeneity, we can observe clear clustering phenomena. However, this should be biased because inhomogeneity make more sense.

Ripley’s K-function

homogeneous

#Estimate the empirical k-function
k_wolf <- Kest(wolf_ppp,correction=c("isotropic"))

k_wolf

## Function value object (class 'fv')
## for the function r -> K(r)
## ................................................................
##      Math.label     Description                                 
## r    r              distance argument r                         
## theo K[pois](r)     theoretical Poisson K(r)                    
## iso  hat(K)[iso](r) Ripley isotropic correction estimate of K(r)
## ................................................................
## Default plot formula:  .~r
## where "." stands for 'iso', 'theo'
## Recommended range of argument r: [0, 341660]
## Available range of argument r: [0, 341660]

plot(k_wolf,
     main = "",
     lwd = 2)

Bootstrapping with 95% CI

E_wolf <- envelope(wolf_ppp,
                  Kest,
                  correction="isotropic",
                  rank = 1,
                  nsim = 19,
                  fix.n = T)

## Generating 19 simulations of CSR with fixed number of points  ...
## 1, 2,  [4:30 remaining] 3,
##  [4:04 remaining] 4,  [3:47 remaining] 5,  [3:31 remaining] 6,
##  [3:21 remaining] 7,  [3:04 remaining] 8,  [2:49 remaining] 9,
##  [2:33 remaining] 10,  [2:19 remaining] 11,  [2:03 remaining] 12,
##  [1:48 remaining] 13,  [1:33 remaining] 14,  [1:17 remaining] 15,
##  [1:01 remaining] 16,  [46 sec remaining] 17,  [30 sec remaining] 18,
##  [15 sec remaining] 
## 19.
## 
## Done.

plot(E_wolf,
     main = "",
     lwd = 2)

For homogeneous proecess, it looks there is clustering.

inhomogeneous

lambda_wolf <- density(wolf_ppp, bw.ppl)
Kinhom_wolf <- Kinhom(wolf_ppp, lambda_wolf,correction="isotropic")
Kinhom_wolf

## Function value object (class 'fv')
## for the function r -> K[inhom](r)
## ................................................................................
##      Math.label              
## r    r                       
## theo {K[inhom]^{pois}}(r)    
## iso  {hat(K)[inhom]^{iso}}(r)
##      Description                                        
## r    distance argument r                                
## theo theoretical Poisson K[inhom](r)                    
## iso  Ripley isotropic correction estimate of K[inhom](r)
## ................................................................................
## Default plot formula:  .~r
## where "." stands for 'iso', 'theo'
## Recommended range of argument r: [0, 341660]
## Available range of argument r: [0, 341660]

# visualise the results
plot(Kinhom_wolf,
     theo ~ r,
     main = "",
     col = "grey70",
     lty = "dashed",
     lwd = 2)

plot(Kinhom_wolf,
     iso ~ r,
     col = c("#046C9A"),
     lwd = 2,
     add = T)

Bootstrapping with 95% CI

lambda_wolf_pos <- density(wolf_ppp,
                          sigma=bw.ppl,
                          positive=TRUE)

#Simulation envelope (with points drawn from the estimated intensity)
E_wolf_inhom <- envelope(wolf_ppp,
                        Kinhom,
                        simulate = expression(rpoispp(lambda_wolf_pos)),
                        correction="border",
                        rank = 1,
                        nsim = 19)

## Generating 19 simulations by evaluating expression  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 
## 19.
## 
## Done.

# visualise the results
par(mfrow = c(1,2))
plot(E_wolf_inhom,
     main = "",
     lwd = 2)

# Zoom in on range where significant deviations appear
plot(E_wolf_inhom,
     xlim = c(0,50000),
     main = "",
     lwd = 2)

After incorporating the assumption of inhomogeneity, significant clustering only appears to exist within around 30000 meters.

Pair Correlation Function

homogeneous

pcf_wolf <- pcf(wolf_ppp)

pcf_wolf

## Function value object (class 'fv')
## for the function r -> g(r)
## ..............................................................
##       Math.label        Description                           
## r     r                 distance argument r                   
## theo  g[Pois](r)        theoretical Poisson g(r)              
## trans hat(g)[Trans](r)  translation-corrected estimate of g(r)
## iso   hat(g)[Ripley](r) isotropic-corrected estimate of g(r)  
## ..............................................................
## Default plot formula:  .~r
## where "." stands for 'iso', 'trans', 'theo'
## Recommended range of argument r: [0, 341660]
## Available range of argument r: [0, 341660]

plot(pcf_wolf)

# visualise the results
plot(pcf_wolf,
     theo ~ r,
     ylim = c(0,20),
     main = "",
     col = "grey70",
     lwd = 2,
     lty = "dashed")

plot(pcf_wolf,
     iso ~ r,
     col = c("#046C9A"),
     lwd = 2,
     add = T)

inhomogeneous

# Estimate g corrected for inhomogeneity
g_inhom <- pcfinhom(wolf_ppp)

# visualise the results
plot(g_inhom,
     theo ~ r,
     ylim = c(0,9),
     main = "",
     col = "grey70",
     lwd = 2,
     lty = "dashed")

plot(g_inhom,
     iso ~ r,
     col = c("#046C9A"),
     lwd = 2,
     add = T)

#Simulation envelope (with points drawn from the estimated intensity)
pcf_wolf_inhom <- envelope(wolf_ppp,
                          pcfinhom,
                          simulate = expression(rpoispp(lambda_wolf_pos)),
                          rank = 1,
                          nsim = 19)

## Generating 19 simulations by evaluating expression  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 
## 19.
## 
## Done.

# visualise the results
par(mfrow = c(1,2))
plot(pcf_wolf_inhom)

# Zoom in on range where significant deviations appear
plot(pcf_wolf_inhom,
     xlim = c(0,60000),
     main = "",
     lwd = 2)

From the pair correlation function, it is evident that clustering is more significant within 10,000 meters. Between 10,000 and 60,000 meters, however, there is an indication of avoidance.

Fit model to predict the spatial distribution of coyotes

Collinearity check

We check the collinearity of variables at first.

cor.im(BC_wolf$Elevation,BC_wolf$Forest,BC_wolf$HFI, BC_wolf$Dist_Water, use = "pairwise.complete.obs")

##             ..1         ..2         ..3         ..4
## ..1  1.00000000 -0.26204718 -0.26625626 -0.03426387
## ..2 -0.26204718  1.00000000  0.06615541  0.04825439
## ..3 -0.26625626  0.06615541  1.00000000  0.13229408
## ..4 -0.03426387  0.04825439  0.13229408  1.00000000

From the correlation matrix above, we can see that the collinearity between Forest, HFI, Dist_Water and Elevation is very weak. All of them can be used for coyote spatial distribution modelling.

Model fitting

Based on the covariates relationships (rhohat plots), we guess that there is an non-linearity relationship between Dist_Water, HFI, Forest and Elevations. The quadratic term of Dist_water and forest might influence of spatial intensity of coyote, and the intensity of coyote might have strong correlation with HFI and Elevations’ higher polynomial terms.

We start at Forest, Dist_Water, and the Dist_Water and Elevation as our start point for ppm building. And utilize AIC as one of our comparsion standard

model_1 <- ppm(wolf_ppp ~ HFI + Forest + Dist_Water + Elevation, data = BC_wolf)

## Warning: Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of
## 2246) of the quadrature points. Occurred while executing: ppm.ppp(Q = wolf_ppp,
## trend = ~HFI + Forest + Dist_Water + Elevation,

model_1

## Nonstationary Poisson process
## Fitted to point pattern dataset 'wolf_ppp'
## 
## Log intensity:  ~HFI + Forest + Dist_Water + Elevation
## 
## Fitted trend coefficients:
##   (Intercept)           HFI        Forest    Dist_Water     Elevation 
## -2.247316e+01  6.903940e+00 -2.058131e-03  9.045978e-06 -1.474970e-03 
## 
##                  Estimate         S.E.       CI95.lo       CI95.hi Ztest
## (Intercept) -2.247316e+01 2.640138e-01 -2.299062e+01 -2.195570e+01   ***
## HFI          6.903940e+00 2.826821e-01  6.349893e+00  7.457986e+00   ***
## Forest      -2.058131e-03 2.492556e-03 -6.943451e-03  2.827188e-03      
## Dist_Water   9.045978e-06 2.622055e-05 -4.234535e-05  6.043731e-05      
## Elevation   -1.474970e-03 1.751079e-04 -1.818175e-03 -1.131765e-03   ***
##                    Zval
## (Intercept) -85.1211503
## HFI          24.4229792
## Forest       -0.8257113
## Dist_Water    0.3449958
## Elevation    -8.4232050
## Problem:
##  Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of 2246) 
## of the quadrature points

AIC(model_1)

## [1] 13741.68

The simple model shows that Forest and Dist_Water might have weak relationship with coyote’s distribution.

model_2 <- ppm(wolf_ppp ~ HFI + Forest + I(Forest^2) + I(Dist_Water) + Elevation, data = BC_wolf)

## Warning: Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of
## 2246) of the quadrature points. Occurred while executing: ppm.ppp(Q = wolf_ppp,
## trend = ~HFI + Forest + I(Forest^2) + I(Dist_Water) +

model_2

## Nonstationary Poisson process
## Fitted to point pattern dataset 'wolf_ppp'
## 
## Log intensity:  ~HFI + Forest + I(Forest^2) + I(Dist_Water) + Elevation
## 
## Fitted trend coefficients:
##   (Intercept)           HFI        Forest   I(Forest^2) I(Dist_Water) 
## -2.240181e+01  6.959465e+00 -1.012476e-02  9.729086e-05 -3.413309e-07 
##     Elevation 
## -1.490738e-03 
## 
##                    Estimate         S.E.       CI95.lo       CI95.hi Ztest
## (Intercept)   -2.240181e+01 2.701828e-01 -2.293135e+01 -2.187226e+01   ***
## HFI            6.959465e+00 2.883950e-01  6.394222e+00  7.524709e+00   ***
## Forest        -1.012476e-02 6.300284e-03 -2.247309e-02  2.223564e-03      
## I(Forest^2)    9.729086e-05 6.939177e-05 -3.871451e-05  2.332962e-04      
## I(Dist_Water) -3.413309e-07 2.736704e-05 -5.397975e-05  5.329709e-05      
## Elevation     -1.490738e-03 1.758718e-04 -1.835441e-03 -1.146036e-03   ***
##                       Zval
## (Intercept)   -82.91350640
## HFI            24.13171099
## Forest         -1.60703315
## I(Forest^2)     1.40205182
## I(Dist_Water)  -0.01247233
## Elevation      -8.47627596
## Problem:
##  Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of 2246) 
## of the quadrature points

AIC(model_2)

## [1] 13741.75

model_3 <- ppm(wolf_ppp ~ HFI + I(Dist_Water) + I(Forest^2) + I(Forest) + Elevation, data = BC_wolf)

## Warning: Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of
## 2246) of the quadrature points. Occurred while executing: ppm.ppp(Q = wolf_ppp,
## trend = ~HFI + I(Dist_Water) + I(Forest^2) +

model_3

## Nonstationary Poisson process
## Fitted to point pattern dataset 'wolf_ppp'
## 
## Log intensity:  ~HFI + I(Dist_Water) + I(Forest^2) + I(Forest) + Elevation
## 
## Fitted trend coefficients:
##   (Intercept)           HFI I(Dist_Water)   I(Forest^2)     I(Forest) 
## -2.240181e+01  6.959465e+00 -3.413309e-07  9.729086e-05 -1.012476e-02 
##     Elevation 
## -1.490738e-03 
## 
##                    Estimate         S.E.       CI95.lo       CI95.hi Ztest
## (Intercept)   -2.240181e+01 2.701828e-01 -2.293135e+01 -2.187226e+01   ***
## HFI            6.959465e+00 2.883950e-01  6.394222e+00  7.524709e+00   ***
## I(Dist_Water) -3.413309e-07 2.736704e-05 -5.397975e-05  5.329709e-05      
## I(Forest^2)    9.729086e-05 6.939177e-05 -3.871451e-05  2.332962e-04      
## I(Forest)     -1.012476e-02 6.300284e-03 -2.247309e-02  2.223564e-03      
## Elevation     -1.490738e-03 1.758718e-04 -1.835441e-03 -1.146036e-03   ***
##                       Zval
## (Intercept)   -82.91350640
## HFI            24.13171099
## I(Dist_Water)  -0.01247233
## I(Forest^2)     1.40205182
## I(Forest)      -1.60703315
## Elevation      -8.47627596
## Problem:
##  Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of 2246) 
## of the quadrature points

AIC(model_3)

## [1] 13741.75

model_4 <- ppm(wolf_ppp ~ HFI*I(Dist_Water) + Elevation, data = BC_wolf)

## Warning: Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of
## 2246) of the quadrature points. Occurred while executing: ppm.ppp(Q = wolf_ppp,
## trend = ~HFI * I(Dist_Water) + Elevation,

model_4

## Nonstationary Poisson process
## Fitted to point pattern dataset 'wolf_ppp'
## 
## Log intensity:  ~HFI * I(Dist_Water) + Elevation
## 
## Fitted trend coefficients:
##       (Intercept)               HFI     I(Dist_Water)         Elevation 
##     -2.271192e+01      7.225534e+00      7.827459e-05     -1.492488e-03 
## HFI:I(Dist_Water) 
##     -1.232840e-04 
## 
##                        Estimate         S.E.       CI95.lo       CI95.hi Ztest
## (Intercept)       -2.271192e+01 2.305704e-01 -2.316383e+01 -2.226002e+01   ***
## HFI                7.225534e+00 3.048343e-01  6.628070e+00  7.822998e+00   ***
## I(Dist_Water)      7.827459e-05 5.509886e-05 -2.971719e-05  1.862664e-04      
## Elevation         -1.492488e-03 1.764431e-04 -1.838310e-03 -1.146666e-03   ***
## HFI:I(Dist_Water) -1.232840e-04 9.712567e-05 -3.136468e-04  6.707884e-05      
##                         Zval
## (Intercept)       -98.503212
## HFI                23.703154
## I(Dist_Water)       1.420621
## Elevation          -8.458754
## HFI:I(Dist_Water)  -1.269324
## Problem:
##  Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of 2246) 
## of the quadrature points

AIC(model_4)

## [1] 13740.81

From those models above, we can see that all terms contains Dist_Water and Forest (no matter what kind it is) has week effect on coyote’s distribution. We remove them from our model.

model_5 <- ppm(wolf_ppp ~ HFI + Elevation, data = BC_wolf)

## Warning: Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of
## 2246) of the quadrature points. Occurred while executing: ppm.ppp(Q = wolf_ppp,
## trend = ~HFI + Elevation, data = list(list(

model_5

## Nonstationary Poisson process
## Fitted to point pattern dataset 'wolf_ppp'
## 
## Log intensity:  ~HFI + Elevation
## 
## Fitted trend coefficients:
##   (Intercept)           HFI     Elevation 
## -22.579650812   7.012690936  -0.001477315 
## 
##                  Estimate         S.E.       CI95.lo       CI95.hi Ztest
## (Intercept) -22.579650812 0.2064665368 -22.984317788 -22.174983836   ***
## HFI           7.012690936 0.2510511123   6.520639798   7.504742075   ***
## Elevation    -0.001477315 0.0001747784  -0.001819874  -0.001134755   ***
##                    Zval
## (Intercept) -109.362278
## HFI           27.933320
## Elevation     -8.452501
## Problem:
##  Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of 2246) 
## of the quadrature points

AIC(model_5)

## [1] 13738.63

Now all terms are effective in our model. And we want to explore if their polynomial fits the model better. Since the higher-degree polynomial will easily make the model fail to converge, we use gam instead of ppm

model_6 <- ppm(wolf_ppp ~ bs(HFI, knots = c(0.2,0.6), degree = 3, df = 4) + bs(Elevation, degree = 3, df = 3), data = BC_wolf, use.gam = TRUE)

## Warning: Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of
## 2246) of the quadrature points. Occurred while executing: ppm.ppp(Q = wolf_ppp,
## trend = ~bs(HFI, knots = c(0.2, 0.6), degree = 3,

model_6

## Nonstationary Poisson process
## Fitted to point pattern dataset 'wolf_ppp'
## 
## Log intensity:  ~bs(HFI, knots = c(0.2, 0.6), degree = 3, df = 4) + 
## bs(Elevation, degree = 3, df = 3)
## 
## Fitted trend coefficients:
##                                       (Intercept) 
##                                       -23.8740545 
## bs(HFI, knots = c(0.2, 0.6), degree = 3, df = 4)1 
##                                         0.2301344 
## bs(HFI, knots = c(0.2, 0.6), degree = 3, df = 4)2 
##                                         5.1913730 
## bs(HFI, knots = c(0.2, 0.6), degree = 3, df = 4)3 
##                                         6.9626038 
## bs(HFI, knots = c(0.2, 0.6), degree = 3, df = 4)4 
##                                         6.5455994 
## bs(HFI, knots = c(0.2, 0.6), degree = 3, df = 4)5 
##                                         8.4684642 
##                bs(Elevation, degree = 3, df = 3)1 
##                                        -6.0917997 
##                bs(Elevation, degree = 3, df = 3)2 
##                                        14.2016702 
##                bs(Elevation, degree = 3, df = 3)3 
##                                       -64.3914811 
## 
## For standard errors, type coef(summary(x))
## Problem:
##  Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of 2246) 
## of the quadrature points

AIC(model_6)

## [1] 13620.24

anova(model_5, model_6, test = "LRT")

## Warning: Models were re-fitted with use.gam=TRUE

## Warning: Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of
## 2246) of the quadrature points. Occurred while executing: ppm.ppp(Q = wolf_ppp,
## trend = ~HFI + Elevation, data = list(list(

## Warning: Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of
## 2246) of the quadrature points. Occurred while executing: ppm.ppp(Q = wolf_ppp,
## trend = ~bs(HFI, knots = c(0.2, 0.6), degree = 3,

## Analysis of Deviance Table
## 
## Model 1: ~HFI + Elevation     Poisson
## Model 2: ~bs(HFI, knots = c(0.2, 0.6), degree = 3, df = 4) + bs(Elevation, degree = 3, df = 3)    Poisson
##   Npar Df Deviance  Pr(>Chi)    
## 1    3                          
## 2    9  6    130.4 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

par_res_elev <- parres(model_6, "Elevation")

## Warning: Some infinite, NA or NaN increments were removed

par_res_hfi <- parres(model_6, "HFI")

## Warning: Some infinite, NA or NaN increments were removed

## Warning: Values for 1 query point lying outside the pixel image domain were
## estimated by projection to the nearest pixel

par_res_forest <- parres(model_6, "Forest")

## Warning: Some infinite, NA or NaN increments were removed

par_res_water <- parres(model_6, "Dist_Water")

## Warning: Some infinite, NA or NaN increments were removed

par(mfrow = c(1,2))
plot(par_res_elev,
     legend = FALSE,
     lwd = 2,
     main = "",
     xlab = "Elevation (m)")
plot(par_res_hfi,
     legend = FALSE,
     lwd = 2,
     main = "",
     xlab = "HFI")

As the AIC and anova report shown, our new GAM model has a better performance now. But it also looks luck performance on the low elevation, and the HFI’s confidence interval has high variance when HFI near to 0.6 or lower than 0.2. We add two knots to capture more detailed variations in coyote activity intensity with respect to elevation, and we increased the weight proportion for HFI > 0.8 because, according to rhohat, the frequency of coyote sightings only begins to rise sharply after HFI exceeds 0.8.

model_7 <- ppm(wolf_ppp ~ bs(HFI, degree = 4, df = 4) + I(HFI > 0.8) + bs(Elevation, knots = c(300,1400),degree = 4, df = 3), data = BC_wolf, use.gam = TRUE)

## Warning: Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of
## 2246) of the quadrature points. Occurred while executing: ppm.ppp(Q = wolf_ppp,
## trend = ~bs(HFI, degree = 4, df = 4) +

model_7

## Nonstationary Poisson process
## Fitted to point pattern dataset 'wolf_ppp'
## 
## Log intensity:  ~bs(HFI, degree = 4, df = 4) + I(HFI > 0.8) + bs(Elevation, 
## knots = c(300, 1400), degree = 4, df = 3)
## 
## Fitted trend coefficients:
##                                              (Intercept) 
##                                            -2.397222e+01 
##                             bs(HFI, degree = 4, df = 4)1 
##                                             2.502930e+00 
##                             bs(HFI, degree = 4, df = 4)2 
##                                             1.058737e+01 
##                             bs(HFI, degree = 4, df = 4)3 
##                                             3.644627e+00 
##                             bs(HFI, degree = 4, df = 4)4 
##                                             8.586450e+00 
##                                         I(HFI > 0.8)TRUE 
##                                            -5.246966e-02 
## bs(Elevation, knots = c(300, 1400), degree = 4, df = 3)1 
##                                             2.285253e-01 
## bs(Elevation, knots = c(300, 1400), degree = 4, df = 3)2 
##                                            -4.159536e+00 
## bs(Elevation, knots = c(300, 1400), degree = 4, df = 3)3 
##                                             9.936512e+00 
## bs(Elevation, knots = c(300, 1400), degree = 4, df = 3)4 
##                                            -3.594385e+01 
## bs(Elevation, knots = c(300, 1400), degree = 4, df = 3)5 
##                                             7.435627e+01 
## bs(Elevation, knots = c(300, 1400), degree = 4, df = 3)6 
##                                            -1.106199e+03 
## 
## For standard errors, type coef(summary(x))
## Problem:
##  Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of 2246) 
## of the quadrature points

AIC(model_7)

## [1] 13604.38

anova(model_6, model_7, test = "LRT")

## Warning: Models were re-fitted with use.gam=TRUE

## Warning: Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of
## 2246) of the quadrature points. Occurred while executing: ppm.ppp(Q = wolf_ppp,
## trend = ~bs(HFI, knots = c(0.2, 0.6), degree = 3,

## Warning: Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of
## 2246) of the quadrature points. Occurred while executing: ppm.ppp(Q = wolf_ppp,
## trend = ~bs(HFI, degree = 4, df = 4) +

## Analysis of Deviance Table
## 
## Model 1: ~bs(HFI, knots = c(0.2, 0.6), degree = 3, df = 4) + bs(Elevation, degree = 3, df = 3)    Poisson
## Model 2: ~bs(HFI, degree = 4, df = 4) + I(HFI > 0.8) + bs(Elevation, knots = c(300, 1400), degree = 4, df = 3)    Poisson
##   Npar Df Deviance Pr(>Chi)    
## 1    9                         
## 2   12  3   21.856 6.99e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

par_res_elev <- parres(model_7, "Elevation")

## Warning: Some infinite, NA or NaN increments were removed

par_res_hfi <- parres(model_7, "HFI")

## Warning: Some infinite, NA or NaN increments were removed

## Warning: Values for 1 query point lying outside the pixel image domain were
## estimated by projection to the nearest pixel

par_res_forest <- parres(model_7, "Forest")

## Warning: Some infinite, NA or NaN increments were removed

par_res_water <- parres(model_7, "Dist_Water")

## Warning: Some infinite, NA or NaN increments were removed

par(mfrow = c(1,2))
plot(par_res_elev,
     legend = FALSE,
     lwd = 2,
     main = "",
     xlab = "Elevation (m)")
plot(par_res_hfi,
     legend = FALSE,
     lwd = 2,
     main = "",
     xlab = "HFI")

Now We can see that the whole plot fitting has improved.

Since in the coyote intensity distribution, we can see that the coyote clustered in the HFI > 0.8. So we see the fitting in the small window.

par(mfrow = c(1,2))
plot(par_res_elev,
     legend = FALSE,
     lwd = 2,
     main = "",
     xlab = "Elevation (m)", ylim = c(-20,20))
plot(par_res_hfi,
     legend = FALSE,
     lwd = 2,
     main = "",
     xlab = "HFI", xlim = c(0.8,1))

We can see that the model fitting well now.

plot(model_7,log=TRUE,se = FALSE,superimpose = FALSE, n= 100)

Here is our residual plot

res <- residuals(model_7)

## Warning: Some infinite, NA or NaN increments were removed

na_indexes <- which(is.na(res$val))
# If you need to exclude NA values explicitly
if (length(na_indexes) > 0) {
  res <- res[-na_indexes]
  plot(res, cols='transparent')
} else {
  plot(res, cols='transparent')
}

From the graph above, we can see that our model can fit the coyote distribution well. # lucking variable check

lurking(model_7, BC_wolf$Forest)

## Warning in LurkEngine(object = object, type = type, cumulative = cumulative, :
## 19 out of 2246 quadrature points discarded because they lie outside the domain
## of the covariate image

## Warning: model was fitted by gam(); asymptotic variance calculation ignores
## this

lurking(model_7, BC_wolf$Dist_Water)

## Warning: model was fitted by gam(); asymptotic variance calculation ignores
## this

There is less lucking variable effect in our model. So it means our variable selected is suitable.

diagnose.ppm(model_7)

## Warning: Some infinite, NA or NaN increments were removed

## Warning: model was fitted by gam(); asymptotic variance calculation ignores
## this

## Warning: model was fitted by gam(); asymptotic variance calculation ignores
## this

## Model diagnostics (raw residuals)
## Diagnostics available:
##  four-panel plot
##  mark plot 
##  smoothed residual field
##  x cumulative residuals
##  y cumulative residuals
##  sum of all residuals
## sum of raw residuals in entire window = -5.324e-10
## area of entire window = 9.483e+11
## quadrature area = 9.388e+11
## range of smoothed field =  [-1.651e-10, 3.034e-10]

We have found an interesting issue that our residuals increase on the southwest coast, where we also have a higher record of coyote observations. Diagnostically, we can see that the peaks of the errors almost all occur in the urban areas of British Columbia: Vancouver, Kelowna, etc. We believe that these prediction errors stem from data sampling. Most of the coyote data relies on human sighting reports, which causes areas with high Human Footprint Index (HFI) to have more coyote records, while areas with less human activity have fewer records, hence our data exhibits sampling bias. Additionally, due to factors like urban environments and mountains, the distribution of coyotes is not homogeneous. This means that our four variables: forest, distance to water, elevation, and HFI, cannot fully encapsulate the factors influencing coyote distribution. To accurately predict coyote appearances, we need more data.